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Abstract 

The aim of this work is to study the numerical solution of the non- 
linear Schrodinger problem using a combination between Witt basis 
and finite difference approximations. We construct a discrete funda- 
mental solution for the non-stationary Schrodinger operator and we 
show the convergence of the numerical scheme. Numerical examples 
are given at the end of the paper 
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1 Introduction 

In this paper we make use of Clifford analysis tools in order to treat a well- 
known partial differential equation of mathematical physics. This treatment 
is based on the work developed by K. Giirlebeck and W. Sprofiig in [T] and 
it is (partially) based on an orthogonal decomposition of the underlying 
function space in terms of the subspace of null-solutions of the correspond- 
ing Dirac operator. While the orthogonal decomposition of Giirlebeck and 
Sprofiig has been applied with success to PDE's such as Lame equations, 
Maxwell equations and Navier-Stokes equations [2J, it works for the station- 
ary case only. 

In [3] an alternative approach was proposed, based on an adding of extra 
basis elements, namely, of a Witt basis. This approach allows the application 
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of the already existent techniques of elliptic function theory developed in 
PQ to time-varying domains. A suitable orthogonal decomposition for the 
underlying function space is, therefore, obtained in terms of the kernel of 
the positive parabolic Dirac operator and its range after application to a 
Sobolev space with zero boundary values. 

After some basic notions about Clifford algebras presented in the next 
section, we will define, in Section 2, a generalization of the parabolic Dirac 
operator introduced in [3] and a generalization of the Teodorescu and Cauchy- 
Bitsadze operators presented in pQ. Moreover, using the previous definitions 
we will obtain a factorization of our equation in terms of basis elements of 
Witt basis and we obtain the fundamental solution for our generic parabolic 
Dirac operator. 

However, the integral representation formulae obtained via this theoret- 
ical method are not suitable for an explicit computation of the solution, 
due to unacceptable convergence rates of the integrals' numerical approxi- 
mation (see pQ for more details). Hence, to avoid this backdraw it becomes 
necessary to study the discrete analogues of the operators, namely discrete 
counterparts for the single- and double-layer potentials. Contrary to differ- 
ence potentials introduced by Ryabenkij [4], where the difference potentials 
are constructed by means of discrete Green functions, we will introduce in 
Section 3 difference potentials based on the discrete fundamental solution. 
An advantage of this approach is that, contrary to discrete Green functions, 
we will obtain an explicit expression for our discrete fundamental solution 
Eh-ir which is independent of the choice or shape of the domain. In Section 
4 we prove the convergence of the discrete counterparts of the analytic op- 
erators introduced in Section 2. This will allow us to establish a convergent 
numerical scheme for the linear non-stationary Schrodinger equation. 

In Section 5 we will adapt the previous algorithm in order to solve nu- 
merically the cubic Schrodinger equation and we will present in Section 6 
some simple numerical examples to show the consistency and stability of our 
algorithm for different mesh sizes h and r. 

2 Preliminaries 
2.1 Clifford algebras 

Consider the n-dimensional vector space W 1 endowed with a standard or- 
thonormal basis {ei, • • • , e n } and satisfying the multiplication rules ej&j + 
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We define the universal Clifford algebra C^o,n as the 2 n -dimensional as- 
sociative algebra with basis given by eo = 1 and = • • • e^, , where 
A = {hi,...,h k } C N = {1, . . . ,n}, for 1 < hi < ■ ■ ■ < h k < n. Each 
element x E C^o,n will be represented by x = Y^a x a^a and each non-zero 
vector x = Yl 1 j=i x j e j e ^ n has a multiplicative inverse given by r^p. We 
denote by x ce °' n the (Clifford) conjugate of the element x £ C£o,n, where 

We introduce the complexified Clifford algebra Ct n as the tensorial pro- 
duct 

C ® Ci^ n = {w = J2 Z ^A, Z A E C, A C iV}, 
A 

where the imaginary unit interact with the basis elements 
l,...,n. 

The conjugation is defined as w = ^2^'ZA C ^A Cio ' n ■ 

We consider the Dirac operator D = Y^j=i e j ~Bx~ which has the property 
of factorizing the n-dimensional Laplacian, that is, D 2 = —A. A C£ n - 
valued function on an open domain Q, u : Q C M n *— > C£ n , is said to be 
left-monogenic if it satisfies Du = on Q. 

Let now Q C M n x M + denote a bounded domain with a sufficiently 
smooth boundary T = dtl, while (0, T), with T > 0, represents its projection 
on the time-domain. A function u : O i— > C£ n has a representation u = 
Y^,A u A e A with C-valued components -u^. Properties such as continuity will 
be understood component-wisely. In the following we will use the short 
notation L p (0,), C k (tl), etc., instead of L p (£l,Cl n ), C k (Q,C£ n ). For more 
details, see [5]. 

Taking into account [3] we will imbed R n into R n+2 . For that purpose 
we add two new basis elements f and f ' satisfying 

f 2 = ft 2 = 0, fft + ftf = 1, fa + e rf = ft ei + erf = 0, j = 1, • • • , n. (1) 

The set {f, ft} is said to be a Witt basis for M. 2 and it will allows us to 
create a suitable factorization of the Schrodinger operator where only partial 
derivatives are used. 

2.2 Factorization of time-evolution operators 

In this section we present a new method for factorizing the Schrodinger 
equation, 

(±id t - A)u(x, t) = 0, (x, t) G n, (2) 
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where C W 1 x R + denotes a bounded domain. For this we will follow the 
ideas presented in [3], [6] and [7]. 

Definition 2.1. For a function u G C 1 (r2) we define the forward (resp. 
backward) parabolic Dirac operator 

D Xj±it u = (D + fd t ± ift) Uj (3) 

where D stands for the (spatial) Dirac operator. 

These operators factorize the correspondent backward/forward time- 
evolution operator ([2]), that is 

{D x , ±u fu = (±id t - A)u. (4) 

We consider now the generic Stokes' Theorem. 

Theorem 2.2. For each u, v G W^(Q), 1 < p < oo it holds 

vda x<t u = \ l(vD x - it )u + v(D Xj+it u)]dxdt, 
an Jo, 

where dcr X) t = {D x + ]dt)\dxdt stands for the contraction of the homogeneous 
operator associated to D x _n with the volume element. 

For the proof of this theorem we refer to [3] . 

We shall construct a fundamental solution for the backward parabolic 
Dirac operator D x _a in terms of a fundamental solution of the backward 
Schrodinger operator. We recall that the function 

/ \ ■ H (t) ( M 2 \ 

e-(x,t) = i- — tt exp j- — (5) 

is a fundamental solution for the backward Schrodinger operator since it 
satisfies 

{-id t - A)e-(x,t) = e-(x,t)(-id t - A) = S(x,t) 
in distributional sense. Therefore, we have 

Definition 2.3. Given a fundamental solution e_ = e_(x,t) for the back- 
ward Schrddinger operator we have as a fundamental solution E- = E-(x,t) 
for the backward parabolic Dirac operator D x _n the function 



E-{x,t) = e-(x,t)D. 



x,—it 



H(t) explitet) / , n x \2 



+ f Ur-*^r +f T • (6) 



{AiTit) n / 2 \ 2t 'V 4 * 2 2t 
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Using the fundamental solution © and the generic Borel-Pompeiu for- 
mula we construct the adequate Teodorescu and Cauchy-Bitsadze operators. 

Definition 2.4. For a function u G L p (Vt), 1 < p < oo, we define the 
correspondent Teodorescu and Cauchy-Bitsadze operators, respectively, as 

Tu(xo,to) = / E-(x — xq, t — to)u(x, t)dxdt, 
Jn 

Fu(x ,t ) = / E_(x - x ,t - t Q )da X)t u(x,t). 
Jdn 

We also have the following decomposition (c.f. [7j). 
Theorem 2.5. The space L p (fl), 1 < p < 2 allows the direct decomposition 

Lp(fi) = Lp(fi) n ker{D xMt ) Z^+it fwjj (fi)^ , 

o 

where W p (O) denotes the space of all functions in the Sobolev space W p (fi) 
uni/i zero-boundary values. 

The previous decomposition of the L p -space allows us to establish two 
projections operators. 

Definition 2.6. Lei 1 < p < 2. VFe define the projectors 

P:L p (n)^L p (n)nker(D X)+tt ) (7) 

and 

Q : L p (n) -> D xMt (w p l (n)^j . (8) 

Theorem 2.7. Let f € L p (f2), /or 1 < p < 2. TTie solution of the forward 
linear Schrddinger problem 

iff - An = f inn 
u = on d£l 

is then given by u = TQTf. 

The proof of this theorem was made in [7] for the case of p = 2. However, 
we remark that it can easily be extended to 1 < p < 2. Moreover, 

1) we can obtain dual results for the backward Schrodinger problem by 
considering a fundamental solution for the forward parabolic Dirac operator 
D Xj+ it on Theorem I2.2( 

2) the above construction can easily be generalized for arbitrary opera- 
tors of type adt — A, where a is a non-zero complex parameter. Indeed, the 
case of a = 1 gives the well-known heat equation while for a = i we have 
the non-stationary Schrodinger equation. 
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3 Discrete fundamental solution for the time-evo- 
lution problem 

3.1 Quaternionic matrix representation of the Witt Basis 

We use the matrix representation of the generators of the real quaternions 
as defined in p], 
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as representatives of a discrete version of the spatial basis for the quater- 
nionic case. 



3.2 Finite differences and time evolution operators 

As already stated we want to investigate a finite difference scheme based on 
the notion of a discrete fundamental solution as described in [8] . We denote 
by 



I fern = (hrm , hm?,, hm^),mi £ Z) and 



{kr, k G Z + } 



equidistant lattices corresponding to space and time discretization, respec- 



tively. For a discrete function u 



g5 H, u(hm, kr) 



u , u 



we have the finite difference approximation for the stationary 



Dirac operators given by 
( 



A 1 " 1 - 



d^U 2 



- d h 3u3 



dl 



-dlu 1 



dlu - d~ 3 u 2 + d h 



dlu 3 



2„,0 



d~ h u 



°-d 
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V 

-v 



( -d, 



"V- 



d~ h 3 u 3 



d~ 2 



V 9, 



h " 

-3„.0 



d^u 1 
dlu' 



dlu 3 
dlu 2 



uDt 



-dlu 1 - d\u 2 
^ V - 



u 

2„,1 



d~ 2 u 3 
dZ\ 3 
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where 

flfV = ^(hrn±he s ,kr)-^(hrn,kr))^ . = ^ ^ ^ j = ^ ^ 3^ 

represent the spatial forward/backward difference operators. We remark 
that these difference Dirac operators factorize the star discretization of the 
Laplace operator, in the sense that 

DtD h + = D-+D+- = -A h e = (lX*<*J e °' 

Moreover, we also have the following (forward) time difference operator 
(see P, 0) 

• u>{hm,(k + l)T) -vP(hm,kT) ■ n Q 

o T u J (hm, kt) = ,,7 = 0, •••,3. 

r 

With the previous definitions we aim to construct a finite difference ap- 
proximation for the parabolic Dirac operators. For this purpose we introduce 
the matrix representations 

f \ ( d T e \ + ( ie \ , , 

where 7+, 7" denote elements which satisfy the following matricial opera- 
tions 

7± ( C D ) = ( -C D ) 7± ' 

(7 ± ) 2 = 0, (10) 
7 + 7~ + 7~7 + = id. 

Using the properties of the previous operators and taking account the 
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multiplication rules (|10p we obtain the following relation 
D 



(D. 



h,±ir) 











d T e 
<9 T e 



7 + ± 



d T e 
d T e 



7 



ieo 
ieo 

ieo 
ieo 
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<9 T e 
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7 T + 
j- + 

7~ + 

d T D- h 




<9 r e 
<9 T e 
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ieo 
ieo 



D 



7 







± 



± 



7 



+ 



7~± 





id T Bo 




9 T e 
9 T e 



o 



id T eo 



= (-A h ± l 9 T ) 



e 
e 



7 
7 + ± 
( 7 + 7 - +7^7+) 

(11) 



i.e., these operators factorize the difference discretization of our time evolu- 
tion operator ([2]). Moreover, due to the fact that the above finite difference 
operators D^ + , Dt~ and d T are approximations of the Dirac operator D and 
of the time partial derivative operator dt, respectively (see [9]), we have that 
([9]) are a finite difference approximations for the parabolic Dirac operators 
D x ±u- 



3.3 Discrete fundamental solutions 

Based on the ideas presented in [9] we introduce the discrete fundamental 
solution for the Schrodinger difference operator —id T — A^ as 



e h - iT (hrn, kr) = iH(kr) (1 + irA h 

where H denotes the Heaviside function and 





5 h (hm), 



(12) 



S h (hm) 



if hm 
if hm ^ 



5 T (kr) 



are the discrete analogues of the Dirac delta function in 
tively. Easy calculations show that, indeed, we have 



i if kr = 
if kr^O 

and M T 



respec- 



(-id T - A h )e h _ iT (hm, kr) = e h _ iT {-id T - A h )(hm, kr) = 5 T (kT)5 h (hm). 

(13) 



Numerical Clifford Analysis for Nonlinear Schrddinger Problem 



9 



By the factorization property (jlip . we have for the discrete fundamental 
solution of the operator D^—ir the function 

Efi, — iT ^h,—ir-^h,—ir- 

Moreover, straightforward calculations give the following matrix repre- 
sentation for the discrete fundamental solution Eh,~i T 

E h _ iT (hm, kr) = 
D^ + eh,-i T \ , a ( e ^ _ + ( e 



D--e h ,- ir ; + eo j7 + - ^ Q ^ j j 

However, it remains to prove that the discrete fundamental solution 
eh -it is indeed an approximation of the fundamental solution ([5]). This 
will be done in the next section. 



4 Discrete operator calculus 

We define the discrete /^-spaces, 1 < p < oo, in the usual way 

g G l p (R 3 h x R+) 

iff 



lfi'll/ p (R3 x R+) 



^ h 3 T\g(hm, rk)\ p \ < oo. 

y(hm,Tk)eM.l xK+ / 

Henceforward, no distinction will be made between the function u : f2 — > 
C 4 and its restriction u = u(hm,kT) to the lattice Jl/i )T = On (M| x K+), 
this distinction being clear from the context. 

4.1 Behavior of the discrete fundamental solution 

We now study the behavior of the discrete fundamental solution (|12|) when 
h and r tend to zero and we prove that it converges in Zi-sense to the 
restriction to the grid of the fundamental solution ([5]). 

Theorem 4.1. Let -r$ < Then for any bounded domain G C M 3 it 

holds 

\\Zh-iT - e—\\h{G h x[Q,+oo) T ) ~^ 

as h,r — > 0. 
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The proof of this theorem is based on |10j . Theorem 1, after adaptation 
to space dimension n = 3 and taking in account that our solutions differ 
from the ones in the case of the heat operator by the relations 



e/i,-ir(") ") = ^ e /i,r(*i *•)> (discrete case). 

Moreover, due to the fact that the constructed discrete fundamental solution 

eh -ir has a conical support domain we obtain the mesh-size condition •™ < 

l' 

67r' 2 ■ 

We remark that Theorem 14.11 implies the /^-convergence of (|12p to ([5]) . 
Also, as an immediate consequence we have 

Corollary 4.2. Under the conditions of Theorem \4-l\ it holds 



for any bounded discrete domain Gh C M 3 , as h,r — > 0. 

While we can prove the convergence of the discrete solution Eh-i T to 
the proofs do not yield the order of convergence due to the nature of 
the continuous fundamental solution of the Schrodinger equation. This will 
be the subject of future work. 

Hence, we can establish the discrete analogues of the Teodorescu opera- 
tor. 

Theorem 4.3. For all u £ lp(^lhr)i 1 < P < +oo, such that u : 0/j T — > C 4 
we have the discrete Teodorescu operator T^-ir satisfying to 




ie(-,i-) (continuous case) 



and 



\\Eh-ir ~ -E , -lb 1 (G h x[o,+oo) T ) 




(14) 



where 






(hn,sr)£ £lh,r 



for all (hm, kr) G ^h,r- 



Proof. We have for T^-ir that 



D. 






(hn,sr)£ £lh,T 
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Since Eh_ iT = eh-i T Dh,-i T and eh -it is a scalar solution, we have 

Dh-irTh-iTu(hm, kr) 
= ^ h 3 T[e h - iT (D h - iT ) 2 (hrn - hn, kr - sr)]u(hn, st) 

(hn,sr)£ £lh,T 

= [$h(hm — hn)5 T (kT — ST)u(hn, st)] 

= u(hm, kr). 

Now we are able to present the following norm estimate. 



□ 



Theorem 4.4. For all u <G l p (^h,r), 1 < P < +oo, such that u : Qh,r — * C 4 
there exists a positive constant C > such that 

\\ T h,-iTu\\i p (n hiT ) < c \\ u \\i p (n htT y 
Moreover, T/j „j r is a continuous operator. 
Proof. Initially we have 

\\ T h,-iTu\\i p (n htT ) = 



r/i 3 \Eh_ iT {hm — hn, kr — sr)u(hn, st)\ p 



< r ^- 3 \Eh,-i T {hrn — hn, kr — st)\ p \u{hn, st) 

\{hn,sT)(z£lh,T 

Let us take C(m, k) = max(/ in sr ) e ^ h r \Eh.—ir(hm — hn, kr — sr)\. Then there 
exists C = max C(m, k) > 0, this maximum being taken over all (m, k) such 
that (hm,kT) € £lh,r, and the result holds. □ 

As we have done for the analytic case we can establish a decomposition 
of the /p-space. 
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Theorem 4.5. For the space / p (f2/j T ), 1 < p < oo, the following direct 
decomposition 

o 

lp{fth,r) = \sex D h _ iT (iatQ hiT ) ® D h - iT (wl(Qh,r)) 
is valid, with correspondent discrete projection operators 
Ph,r ■ lp{fth,r) ^->kerD h _ iT (mtQ h , T ), 

O 

Qh,r ■ l P (^h,r) >-> D h _ iT (vjp(Q h , T )), 

o o 

where Wp(Qh,r) denotes the discrete counterpart of the Sobolev space WJ;(Q). 

4.2 Convergence of the discrete operators 

We say that u € C l,a (Q) if its first derivatives are a-H61der continuous. 

Theorem 4.6. Let u £ C 1,a (£l). Then it holds T^_^ T u — ► Tu as h,T tend 
to zero. 

Proof. In order to prove the above result we introduce the regularized Teodor- 
escu operator (see [TT] ) 



T £ u(x,t) = / E £ _(x — z,t — r)u(z,r)dzdr, 
Jn 

where 

E e _{x,t) = e~ € ^E_(x,t) 

stands for a regularization of the fundamental (continuous) solution E_ 
and, therefore, it converges in the sense of tempered distributions to £7_ as 
e —* 0. In a similar way, we construct the regularized discrete operator Tjj _ iT 
in terms of the discrete analogue of the regularized fundamental solution 

I hm,]^ 

E £ K _ iT = e" £ — fi Ar « T . 
By definition, we have 

\Th _ ir u(hm, kr) — T £ u(hm, kr)\ 



< 



E £ h _ ir {hm — hn, kr — sr)u(hn, sr)h?T 

{hn,sT)&lh, T 

— I E £ _(hm — z,kr — r)u(z, r)dzdr 
Jn 



(16) 
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Due to the singularity of the continuous fundamental solution £71, we 
will split the continuous domain f2 into parallelepiped W(hn, st) centered at 
the points (hn, st) of the lattice Vth T with side-lengths h and r, respectively. 
Furthermore, let p, q £ N be such that | + | = 1. We have then 

(del) < 



J _ ir (hm — hn, kr — st) — E e _(hm — hn, kr — sr)]u(hn, sr)/i 3 r| 



+ 



\E £ _{hm — hn, kr — sr)u(hn, sr)h 3 T 

(fen, sr) Efi^ T 

- / E £ _{hm — z,kr — r)u(z, r)dzdr } 

Jw(hn,sr) 



(17) 



We use Holder's inequality on the first term and by a convenient adding up 
we get 



(HID < \\E h _ iT - E_\\t p (n hiT) \\u\\ lq{Qhr) 
+ I E £ _(hm — hn, kr — sT)[u(hn, st) — u(z, r)]h 3 r\ 

(hn,sr)S:fih n zGW (hn,sr) 



+ I | \E e _ (hm — hn, kr — st) — E e _ (hm — z, kr — r)] u(z, r) | dzdr ] . 

IW(hn,ST) 



J 

Jw 



(l2(hn,sr)) 



For the term (Ii) we obtain 



(h) < \E £ _(hm - hn, kr - st) 

x / \u(hn, st) — u(z, r)\dzdr 

Jw(hn,sr) 

< \Et(hm — hn,kT — st)\ 

(hn,ST)Gf2|j T , zdW(hn,sr) 

xC / \(hn — z,ST — r)\ a dzdr, 

Jw{hn,sr) 



which goes to zero as h, r — > 0. 
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Finally the term (l2(hn, st)) can be estimate using its Taylor series ex- 
pansion and Holder's inequality 

(h{hn,sT))< / | [Et{hm — hn, kr — st) — E £ _(hm — z, kr — r)] u(z,r)| 

JW(hn,ST) 

< / \\7 E £ __(hm — z,kT — r) ■ (hn — z, st — r) \u(z, r)\dzdr 

Jw(hn,ST) 

< \\VEt(hm- -,kT - ■) ■ (hn--,ST - ^\\L q {W(hn,ST))\\u\\ Lp {W(hn,ST))i 

and again we have that Y,(hn,sr)en h T , zeW(hn,sr)ih{hn, st)) goes to zero as 
h,T 0. 

Hence, by e — > we obtain convergence of the discrete Teodorescu op- 
erator Th-ir to the continuous one. □ 

Moreover, we notice that we have convergence in l p , 1 < p < oo, of the 
regularized discrete Teodorescu operator _ iT to the regularized continuous 
operator T e . 

We now prove the convergence of the discrete Cauchy-Bitsadze operator 
F^-ir = I — Th-irDh^T-. Moreover, in what follows we will consider the 
sub-domains ft* = {x G R 3 : (x, t) G ft} and fP = {i G K+ : (x, t) G ft} . 

Theorem 4.7. // u G kerD^.^ zs suc/i £/ia£ u G C 1,0 (ft) /or some < a < 
1 then we have 

\\u - F h _ iT u\\i p(nhT) < C\\u\\ c i, a{n) (h a + T a ), 

for a positive constant C > 0. 

Proof. We use the definition of Fh_i T , Theorem 14.41 and the fact that u G 
ker .D^-jt. We get then 

\\u - Fh-i T u\\i p (n h ^ = \\T h -i T D h - iT u\\i p (n KT ) 

= \\Th-ir(F) h -i T U - D x _ it u)\\i p{ Q hT ) 

< Ci | \D h - iT u - D x _ it u\ \i p (n h T ) 

< d (\\D h u - D x u\\ lp{Qh T) + \ \d T u - d t u\\ lp{Qh T ^ 



< c x 



\Dh.u(hm, kr) — D x u(hm, kr)\ p h 3 T 

(hm,kr)£Qh t 



+ [ \d T u(hm,kT)-d t u(hm,kT)\ p h 3 T\ p., 

(hm,kr)£Qh T 
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Additionally, we remark that u G C l,a {Q) implies both 

U(;t) G C 1,a (n*), u(x, •) G C 1 '"^). 

Moreover, we have (c.f. [I], p. 268) that 

\D+'u(hm,kr) - D x u(hm,kT)\ < K{kr)\\u{-, kT)\\ c i, a{nkT) h a , (19) 

a similar result holding for D^ + , and 

\d T (hm,kT) — dtu(hm,kT)\ < K(hm)\\u(hm,-)\\ c i, a ^Qhrn-jT a , (20) 

for some positive constants K(kr), K(hm). Using these two inequalities we 
have 



HD < Ci 



^ l/p 

£ Kf(fcr)|K.,fcr)||^ 0(nfcr) W"fc s r l 

(/im,fcT)Gf2h lT 

'/•' 

+ I Y, ^(hm)\\u(hm r )\\ P cha{nh ^rP a h^ 

We now take iT = maxj) h r {.ftr(fcr), K(hm)} > and we recall that 
\\u(hm,-)\\ c i, a( Qhm.) < \M\ch<*(si), \\u(; kT)\\ c i, a{n k T) < \\u\\ c i,*( U y 
Hence 

(HED < CxKVoi{n h , T )\\u\\ c ^ m {h a + r a ). 

□ 

We are now in conditions to prove the convergence of the discrete pro- 
jection operator Qh, T to its continuous counterpart ([8]). 

Theorem 4.8. Let u G L P {Q) for some 1 < p < oo. Then it holds for the 
projector Qh, T 

\\Qh,rU - Qu\\i p{ Q hr) -> as h,r -^0 
for a positive constant C. 
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Proof. We start from the equality 

Qh,rU~Qu = Qh,r(Pu + Qu) - Q(Pu + Qu) 
= Qh,rPu + Qh,rQu ~ QU 

and we wish to obtain estimates for the terms Qh, T Pu and (Qh,r — I)Qu (we 
recall that, being projection operators, Q(Pu) = and Q 2 = Q). 

Since Pu = FPu and Qh, T Fh,-i T u = 0, for the first term we obtain 

Qh, T Pu = Q hiT FPu - Qh^Fh-irPu 

= QhA F ~ Ph,-ir)Pu 

= Qh,r(I ~ Fh-ir -TD x - it )Pu 

= QhA 1 - Fh-ir)Pu 

and, therefore, by Theorem 14.71 we get the following estimate 

\\Qh,TPu\\l p (n h>T ) < \\Qh,r\\ \\ Pu ~ F h-irPu\\l p {Q h T ) 

< C\\Q h , T \\ \\Pu\\ c i, am (h a + r a ), 

taking in account that Qh, T has bounded norm. Moreover, due to the fact 
that P is the projection into the kernel of Dh-i T , it holds ||-P'u||c ,1 > a (Q) < °°- 
For the second term we remember that Qu can be written as Qu = 

o 

D x _ng where g GM^ 1 (&■)■ This leads to 
(Qh,r ~ P)Qu = (Q h ,r ~ I)D x ^ it g 

= Qh,r(Dx,-itg - F>h-ir9) + Qh,rP>h-ir9 ~ P>x~U9 

= QhA D x~itg - D h _ iT g) + (D h ^ ir g - D x _ it g), 

since Qh, T Dh,-irg = Dh,-ir9- Hence, taking into account the previous cal- 
culations, Theorem 14.71 and relations (|19p and (|20p we finally obtain 

\\(Qh,r ~ P)Qu\\i p (n h>T ) < {\\QhA\ + l )\\ D h-ir9 ~ Dx,~it9\\i p (n hiT ) 
as h, t goes to zero. □ 

The above discrete operators allow us to establish a discrete equivalent 
of Theorem 12.71 

Theorem 4.9. Let f £ h(^h,r)- The solution of the discrete Schrddinger 
problem 

( (id T - &h)u = f in tt hjT 
\ u = on dQh,T 

is given by u = T h _ iT Q Kr T h _ iT f. 
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5 The non-linear Schrodinger problem 

Let us now consider the non-linear Schrodinger problem 

idtu — Au = M(u) in Q 
u = on d£l 

where M(u) = \u\ 2 u + f, with / G L 2 (0), and \u\ 2 = Ej=o( ui ) 2 - This 
problem can be reduced to 

u = TQTM(u) in ft, (21) 

a problem for which the next theorem proves existence and uniqueness of 
solution (see [6], [7] for details). 

Theorem 5.1. The problem h21\) has an unique solution given in terms of 
the iterative method 

u n+ i = TQTM{u n ) 
if f £ L<2 (f2) satisfies the condition 

1 

L 2 - 36 . 2"i+i ' 

o i 

Moreover, the iteration method converges for each starting point uq Giy 2 (^) 
such that 

IMk, < 6 , ln+1 + w > 



1 1 f I I 

With W = y' 36 , 2 2(m + l) - ~2rnTT" • 

Based on the discrete operators previously introduced we construct the 
discrete version of problem (I21h for our bounded domain 

u = T h - iT Q hiT T h - iT M(u) in tt h>T . (22) 

Indeed, let v be a solution of ([22]) . Then 

(i9 T - A/j)u = D h ^ iT D h -i T [T h -i T Q h)T T h _ iT M(v)] 
= D h ^ iT [Q h>T T h _ iT M(v)} 
= M(v), 

and due to the properties of the projector Qh, T w e have v = on dQh,r- 

Using the same ideas as in the continuous case (see [7]) we get results 
regarding the convergence and uniqueness of the discrete iterative method 

Un+l = Th-i T Qh,TTh~i T M(u n ). 
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Theorem 5.2. If f € h(^h,r) then the discrete problem |HJ) has a unique 

o 

solution u (Lw\ (^h,r) whenever 

\\f\\h(n h>T ) < 36ChT 

o 

and the initial term u$ £lw\ iP'hr) satisfies 

H«o||i3(n fc , T ) < ^ \-Wh, T , 



The proof of this theorem, being similar to the one in the continuous 
case, will be omitted. 

The following result shows that the solution obtained for the discrete 
problem, which we will denote by u*, converges to the solution obtained for 
the continuous, which we will denote by u. In the proof of the following 
theorem the restriction of M(u) to the space-time grid will be denote by 
M Kt {u). 

Theorem 5.3. Let f € £2(0). Then u* converges to u in VLf lT whenever 
h,r ->0. 

Proof. Again, we need to use the regularized Teodorescu operator. We shall 
denote u% = T£ _ iT Q hjT T^_ iT M h>T (ul) and u e = T e QT e M(u e ). We have 

\K-u% 2{nhT) < \\T^ iT Q hjT T^_ iT M h;T (u e ) -T*QT^M(u^\\ l2{nhT) 

V v — ' 

(I) 

+ \\T^_ iT Q hjT T^_ iT (M htT (ul) - M htT (u e )) \\ l2(nh T) 
< (I) + C KT \\ul - u e \\ l2{ n hT) (\K\\ h( a htT ) + W u % 2 (n h}T ) 

which implies that 

i-^, T (||<|| i2(Chir) + ||^|| i2( ^ T) xl 



- u 



\h(n h , T ) < (I) 



where Ch, T is a positive constant which depends from h and r. By Theorem 
5.2 we can guarantee that 
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This inequality, together with Theorem 15.11 ensures that for sufficiently 
small h and r, the following relation 

1 ~ Ch, T (\K\\ H n h , T ) + \\ u % 2 (n h , T )) > 

holds. Therefore, the convergence of u* to u depends only on the term (I). 
Hereby, we have 

(I) = \\Tl_ lT Q KT Tl_ lT M K M) ~ QT*M(u^\\ h(nhT) 

< \\n,-i T Q h ,rn,-rr K*r(0 ~ lli a (fl fc T ) 

V v ' 

(A) 

+ \\U,-irQh,r {U,-ir ~ T*) M*(« e )| 1^,.) + \\U,. lT (Qh,r ~ Q) T* M* (u*)]^^ 

(B) (C) 

+ \\ T h-irQh,r {Th-ir ~ Tt ) f\\l 2 (n hT ) + \\ T h-ir (Qh,r ~ Q) f\\l 2 (Sl h t )j 

, — ' » , — ' 

(D) (E) 

where M*{u e ) = |u e | 2 u e and M£ T (u e ) denotes its restriction to the space- 
time grid. By Theorem 4.6 we can say that (B) and (D) tend to zero as 
h, t — > 0. Also, Theorem 4.8 implies the same result for both (C) and (E). 
Finally, for (A) we have, from the boundedness of the discrete operators, 
the following relation 

\\n,- lT Qh,rT? h _ iT (M h * T (u e ) - M*(u e )) || Ja( „ fciT) 

< H^-*rOA.r^-irlll a (n fc , r )l|M^ T (ti e ) - M* (u^^ 

< ClCh,r, 

where C\ is a finite constant and Ch, T is a constant which depends on h and 
r and goes to zero with h and r. Therefore, (I) tends to zero when h, r — ► 0, 
thus, proving our result as e — > 0. □ 



6 Numerical Examples 

In order to study the rate of convergence of our method for different mesh 
sizes, we shall present some numerical examples. For simplicity sake, we 
shall use a cubic space domain [—a, a] 3 with an equidistant discretization 
grid of (N + l) 3 points. Also, for the discretization of the time domain we 
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shall consider an equidistant grid with M+l mesh-points. At this point, we 
emphasize that the choice of M and N takes into account the restriction 
< imposed by Theorem 4.1. 

For all the examples below we will be presenting a table with the I 1 — error 
between the approximated solution and the exact solution at given instants 
of time. 

Example 1: As a first example, we consider an exact real- valued C°° 
solution u = (0, tti, U2, U3) for the problem (|2"Tj) . where 

ui(x,t) = e~ Xl cos (jrt + — J sin^TraJia^s) 

u 2 (x,t) = u 3 (x,t) = 0, 

and the corresponding right hand side / = idtu — An — |u 2 |ii. 

In the following table we show the approximation error between the 
exact solution u and its discrete approximation Uh )T on the domain = 
[—5, 5] 3 x [0,2] for different mesh sizes. 



Table 1 



N 


M 


t=0 


t=0.4 


t=0.8 


20 


450 


2.3313x10- 


-3 


1.2799x10 


-3 


5.7386x10- 


-4 


25 


703 


1.5265x10" 


-3 


8.3774x10 


-4 


3.7642x10- 


4 


30 


1013 


1.0765x10" 


-3 


5.9073x10 


-4 


2.6569x10- 


4 


35 


1378 


7.9982x10- 


-4 


4.3844x10 


-4 


1.9706x10" 


4 


40 


1800 


6.1732x10" 


-4 


3.3895x10 


-4 


1.5228x10" 


4 


45 


2278 


4.9075x10- 


-4 


2.6919x10 


-4 


1.2107x10" 


4 


50 


2813 


3.9937x10- 


-4 


2.1923x10 


-4 


9.8534x10- 


-5 


55 


3404 


3.3132x10- 


-4 


1.8193x10 


-4 


8.1714x10" 


-5 


N 


M 


t=1.2 


t=1.6 


t=2 


20 


450 


2.5728x10" 


-4 


1.1633x10 


-4 


5.3040x10" 


-5 


25 


703 


1.6914x10- 


-4 


7.5998x10 


-5 


3.4520x10- 


-5 


30 


1013 


1.1950x10- 


-4 


5.3548x10 


-5 


2.4266x10- 


-5 


35 


1378 


8.8572x10- 


-5 


3.9810x10 


-5 


1.7992x10" 


-5 


40 


1800 


6.8416x10- 


-5 


3.0738x10 


-5 


1.3868x10" 


-5 


45 


2278 


5.4362x10- 


-5 


2.4450x10 


-5 


1.1014x10" 


-5 


50 


2813 


4.4226x10- 


-5 


1.9878x10 


-5 


8.9580x10- 


-6 


55 


3404 


3.6700x10- 


-5 


1.6502x10 


-5 


7.4280x10" 


-6 



l\— error between the approximated solution and the exact solution 



at different instants 

The following graphics (Figures 1. and 2.) show the evolution of the 
l\— norm for the approximation error, with respect to the space-mesh and 
to the time-mesh, respectively. 
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x10 J X10" 




20 25 30 35 40 45 50 55 50u 10 Q0 1500 2000 2500 3000 

Number of points of the space grid Number of points of the time grid 



Figure 1: ? x -error for different space Figure 2: /i-crror for different time 
steps. steps. 

Example 2: In this example we consider an exact complex- valued C°° 
solution u = (0,ui,U2,u^) of (|2~T|) . where 

ut{x,t) = (e~* - l) {x\ - 25) {x\ - 25) {x\ - 25) 
u 2 {x,t) = 0, uz(x,t) = (e~* - l) sin(7rxix 2 2;3) e lXlt . 

Below is the table with the error of approximation between the exact 
solution u and its discrete approximation Uh T on the domain Q = [—5, 5] 3 x 
[0,2], for different mesh sizes, 
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Table 2 



AT 

IN 


TV T 

M 


t=0 


t=0.4 


t=0.8 


20 


450 


4.8846x10 


-3 


2.6819x10' 


-3 


1.2024x10" 


-3 


25 


703 


3.1692x10 


-3 


1.7323x10 


-3 


7.8152x10" 


4 


30 


1013 


2.2183x10 


-3 


1.2172x10 


-3 


5.4746x10" 


4 


35 


1378 


1.6404x10 


-3 


8.9923x10 


-4 


4.4166x10" 


4 


40 


1800 


1.2613x10 


-3 


6.9250x10 


-4 


3.1112x10" 


4 


45 


2278 


9.9988x10 


-4 


5.4847x10 


-4 


2.4668x10" 


4 


50 


2813 


8.1181x10 


-4 


4.4563x10 


-4 


2.0029x10" 


-4 


55 


3404 


6.7227x10 


-4 


3.6914x10' 


-4 


1.6580x10" 


-4 


N 


M 


t=1.2 


t=1.6 


t=2 


20 


450 


5.3907x10' 


-4 


2.4374x10 


-4 


1.1140x10" 


-5 


25 


703 


3.5116x10 


-4 


1.5779x10 


-4 


7.1668x10" 


-5 


30 


1013 


2.4623x10 


-4 


1.1033x10 


-4 


5.0000x10" 


-5 


35 


1378 


9.0828x10 


-4 


8.1648x10 


-5 


3.6900x10" 


-5 


40 


1800 


1.8166x10 


-4 


6.2800x10 


-5 


2.8334x10" 


5 


45 


2278 


1.1076x10 


-4 


4.9816x10 


-5 


2.4428x10" 


-5 


50 


2813 


8.9900x10 


-5 


4.0406x10 


-5 


1.8210x10" 


-5 


55 


3404 


7.4468x10 


-5 


3.3484x10 


-5 


1.5072x10" 


-5 



l\— error between the approximated solution and the exact solution 



at different instants 

followed by the graphics (Figures 3. and 4.) of the evolution of the approx- 
imation error for the correspondent space and time mesh sizes considered. 




20 25 30 35 40 45 50 55 50 o 1000 1500 2000 2500 3000 

Number of points of the space grid Number of points of the time grid 



Figure 3: h -error for different space Figure 4: h -error for different time 
steps. steps. 
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Example 3: Finally, we conclude with an example of an exact solution 
of lower regularity on the domain f2 = [—5, 5] 3 x [0,2], namely an exact 
C 1 — solution u = (0,ux,U2,us) of (f21~j) . with 

ux(x, t) = (e~* - 1) (g(xi) - g(-xx)) (g{x 2 ) - g{-x 2 )) (g(x 3 ) - g(-x 3 )) 

u 2 (x,t) = u 3 (x,t) = 0, 



where g is the auxiliary B-spline of order 3 

6 



g(y) 



3 2 ' 2 24 



11 , y _ 7(y~2) 2 3fa-2) 3 
24 " r 8 8 8 



U _ V , fa-3) 2 
6 2 4 
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if < y < 1 

if 1 < y < 2 

if 2 < y < 3 

if 3 < y < 5 



Again, the corresponding right hand side / = idtu — Au — \u 2 \u. The 
following table gives the error of approximation between the exact solution 
u and its discrete approximation n/j r for different mesh sizes considered. 

Table 3 



N 


M 


t=0 


t=0.4 


t=0.8 


20 


450 


5.0846x10" 


-3 


2.8819x10 


-3 


1.4024x10- 


-3 


25 


703 


3.7149x10- 


-3 


2.0388x10 


-3 


9.1607x10- 


-4 


30 


1013 


2.7242x10" 


-3 


1.4948x10 


-3 


6.7232x10- 


-4 


35 


1378 


1.9355x10- 


-3 


1.0610x10 


-3 


4.7688x10- 


-4 


40 


1800 


1.4763x10- 


-3 


8.1058x10 


-4 


3.6402x10- 


-4 


45 


2278 


1.1856x10- 


-3 


6.5030x10 


-4 


2.9248x10- 


-4 


50 


2813 


9.1813x10- 


-4 


4.5629x10" 


-4 


2.0291x10- 


-4 


55 


3404 


8.0086x10- 


-4 


4.3975x10' 


-4 


1.9751x10" 


-4 


N 


M 


t=1.2 


t=1.6 


t=2 


20 


450 


7.3907x10" 


-4 


2.4437x10' 


-4 


1.9111x10" 


-4 


25 


703 


4.1162x10- 


-4 


1.8495x10 


-4 


8.4088x10- 


-5 


30 


1013 


3.0239x10- 


-4 


1.3550x10 


-4 


6.1402x10- 


-5 


35 


1378 


2.1434x10- 


-4 


9.6336x10 


-5 


4.3534x10- 


-5 


40 


1800 


1.6362x10- 


-4 


7.3510x10 


-5 


3.3166x10- 


-5 


45 


2278 


1.3133x10- 


-4 


5.9066x10 


-5 


2.6610x10- 


-5 


50 


2813 


9.9006x10- 


-5 


4.4078x10 


-5 


1.9410x10" 


-5 


55 


3404 


8.8712x10" 


-5 


3.9888x10 


-5 


1.7956x10" 


-5 



Zi— error between the approximated solution and the exact solution 
at different instants 
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The next graphics (Figures 5. and 6.) show the evolution of the approx- 
imation error in l\— norm for the different space mesh size and time mesh 
size considered. 




10 35 40 45 50 55 500 1000 1500 2000 2500 3000 

Number of points of the space grid Number of points of the time grid 



Figure 5: h -error for different space Figure 6: h -error for different time 
steps. steps. 

Taking into account the previous graphics we are able to observe that 
the order of convergence for the space coordinate is, in all the examples, of 
order 0(h 8 ), while for the time coordinate we get, in all the examples, an 

3 

order of convergence of order 0(t^). We remark that our method seems to 
be stable under functions of lower regularity, since the order of convergence 
for the space and time coordinates remains same in all the three examples. 
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